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Abstract. 

We study hexagon patterns in non-Boussinesq convection of a thin rotating layer 
of water. For realistic parameters and boundary conditions we identify various linear 
instabilities of the pattern. We focus on the dynamics arising from an oscillatory 
side-band instability that leads to a spatially disordered chaotic state characterized 
by oscillating (whirling) hexagons. Using triangulation we obtain the distribution 
functions for the number of pentagonal and heptagonal convection cells. In contrast 
to the results found for defect chaos in the complex Ginzburg-Landau equation and 
in inclined-layer convection, the distribution functions can show deviations from a 
squared Poisson distribution that suggest non-trivial correlations between the defects. 
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1. Introduction 

Rayleigh-Benard convection has served as an excellent paradigm to study spatio- 
temporal chaos in pattern-forming systems. For small Prandtl numbers, for which 
large-scale flows driven by the curvature of the roll pattern are important, quite 
complex dynamical roll patterns exhibiting various types of defects like spirals, targets, 
dislocations, and disclinations are found [1]. Another type of spatio-temporally chaotic 
state is obtained if the system is rotated about a vertical axis. Then the Kiippers-Lortz 
instability [2] leads to the formation of domains of rolls of different orientation, which 
are separated by domain walls. Due to the propagation of the domain walls or by the 
nucleation of defects in the bulk of the pattern a persistent switching of the local roll 
wavevector of the different domains is induced [3]. 

While striped planforms typically prevail in convection, hexagonal cellular patterns 
are observed at onset if the temperature variation across the layer is large and the 
variation of the material parameters with the temperature becomes substantial. In this 
case the Boussinesq approximation is invalid and a certain up-down reflection symmetry 
is strongly broken [4]. Recently, results have been presented also for spatio-temporal 
chaos based on a hexagonal rather than a stripe-like planform. Within the framework of 
Ginzburg-Landau equations and of order-parameter equations of the Swift-Hohenberg 
type the effect of rotation [5, 6, 7] and of large-scale flows [8] on hexagonal patterns has 
been investigated. Two interesting scenarios of spatio-temporal chaos have been found 
that are both due to the rotation. One scenario arises from the Hopf bifurcation to 
oscillating hexagons that replaces the common steady transition from hexagons to rolls 
[4] in the presence of rotation [9, 10]. The oscillations of the hexagons are described 
by the complex Ginzburg-Landau equation [11]. If the transition to the oscillating 
hexagons occurs close to threshold the parameters of the complex Ginzburg-Landau 
equation are generically in a regime in which spatially homogeneous and spatially 
modulated oscillations (i.e. traveling waves in the oscillation amplitude) are bistable 
with spiral-dominated defect chaos [6]. For general initial conditions, which induce 
the formation of spirals, the oscillations of the hexagons are therefore expected to 
be spatio-temporally chaotic. The other scenario involves the steady- hexagon state 
below its transition to oscillating hexagons. Within a weakly nonlinear Swift-Hohenberg 
framework a chaotic state was found that is maintained by the nucleation of dislocations 
in the vicinity of already existing penta-hepta defects. Due to this induced nucleation the 
probability distribution function for the number of defects in the pattern is substantially 
broader than that obtained in simulations of defect chaos in the complex Ginzburg- 
Landau equation [12] and in a model for Ca + -waves [13], as well as in experiments on 
electroconvection of nematic liquid crystals [14] and on undulation chaos in inclined- 
layer convection [15, 16]. 

Motivated by the interesting dynamics of defects in the oscillation amplitude or 
in the hexagonal pattern itself that were obtained in the Ginzburg-Landau and the 
Swift-Hohenberg models we investigate here the dynamics of hexagonal patterns in the 
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presence of rotation on the basis of the standard hydro dynamical description of non- 
Boussinesq convection [17, 4, 18, 19]. Thus, we use the Navier-Stokes equations coupled 
to the heat diffusion equation (called NONBE henceforth, see below) corresponding to 
convection in a layer of water near its density maximum. We perform Galerkin stability 
calculations for the hexagonal pattern and direct numerical simulations of the NONBE 
and find spatio-temporal chaos that involves an oscillatory instability of the hexagons as 
well as defects in the hexagonal background pattern. In contrast to the defect chaos of 
the oscillating hexagons studied in the weakly nonlinear regime [5], here the oscillatory 
instability does not correspond to spatially homogeneous oscillations but involves a 
spatial modulation of the oscillation amplitude and it leads to the nucleation of defects 
in the hexagonal background pattern. 



2. Basic Equations 

A horizontal fluid layer of thickness d that is heated from below and cooled from above 

— * 

and that is rotated about a vertical axis with angular velocity cuh is described by the 
Navier-Stokes equation for the momentum, 

dt(pUi) + djipujUi) = - dip- pg5 i3 + 

dj (up (diUj + djUi)) + 2pue Um kiu m , (1) 

the continuity equation, 

d t p + d j ( P u j ) = 0, (2) 
and the heat equation, 

dt T + u 3 d 3 T = ±-d 3 (\d ] T). (3) 

pc p 

Here u = (ui, 112,113) is the fluid velocity, T the temperature, p the density, p the 
pressure, g the acceleration of gravity, v the viscosity, A the heat conductivity, c p the 
specific heat of the fluid, and k is the unit vector in the z-direction. The Kronecker 
delta is given by 5ij and is the unit antisymmetric tensor of the rank 3. We assume 
the Einstein summation convention. We have omitted viscous heating, volume viscosity, 
and the centrifugal force. 

In our analysis we take realistic rigid boundary conditions at the top and bottom 
plates and keep the top and bottom temperatures fixed, 

u = at z = (4) 
T = T + ^f a tz = -± (5) 

T = T - — at z = +-. (6) 

Here T denotes the mean temperature and AT > is the temperature difference 
across the layer. In the lateral directions we chose periodic boundary conditions or 
introduce a suitable ramp in the imposed temperature gradient (see below). 
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In this paper we focus on weakly non-Boussinesq convection. We therefore keep 
the temperature dependence of the various fluid properties to leading order and expand 
them about the mean temperature To in the form [4] 

^=1+73^+ (9) 
A +73 AT + - W 

cJT) T-T 

-^ = l+ 74 — -^ + ... (10) 
c p0 Ai 

where p , u , A , and c p0 denote the values of the respective quantities at the mean 
temperature T . The dots denote higher-order terms to be neglected in the sequel. 
In line with a clever experimental strategy we assume that T is kept fixed when 
changing the main control parameter AT. Thus, the quantities 7j/AT, % = 0,2,3,4, 
which describe, respectively, the slope of the density, viscosity, heat conductivity, and 
heat capacity at T , are fixed as well. This implies that the 7^, which give the change 
of the respective fluid property across the fluid layer, are linear in the main control 
parameter AT and therefore need to be adjusted with increasing AT (see below). The 
usual heat expansion coefficient at T = T is given by a = 70/ AT, but going beyond 
the Boussinesq approximation also the curvature of p(T) at T , which is proportional to 
7071/AT 2 , comes into play. 

Note that for constant fluid parameters the volume viscosity term and the 
centrifugal force can be absorbed into the pressure term. This is not true when the 
properties depend on the temperature. While the ^-dependence of the volume viscosity 
generates only a higher-order correction, neglecting the centrifugal force is only possible 
as long as the rotation rates are not too large for a given system size L, i.e. as long as 
joLu 2 /g is small. 

To render the equations dimensionless we use the usual scales for the length (d), time 
(cP/ko), velocity (kq/oI), pressure (uoHoPo/d 2 ), and temperature {vqKo/ a^gd 3 = AT/R). 
Here we have introduced the heat diffusivity k = Ao/po c />o- This leads to the usual 
dimensionless quantities like the Prandtl number Pr = uq/k,q, the Rayleigh number 
R = a ATgd 3 /uqHq, and the rotation rate Q = ujd 2 /u Q . Finally, we write the equations 
in terms of the dimensionless momenta Vi = pdui/ p n instead of the velocities. 

Since the fluid velocities are small compared to the sound velocity we make the 
anelastic approximation [20] (see also [18]) and neglect the time derivative in the 
continuity equation (2). This simplifies the computation considerably since it reduces 
the number of evolution equations. Furthermore, i>j becomes a solenoidal field, which 
can be represented in the standard poloidal-toroidal decomposition by two velocity 
potentials [18], which automatically enforces the mass conservation. 
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The conduction solution {y = 0) of (3) with (5,6,7,9,10) is given by 

T cond = T q + r(-z- ^-(z 2 - i) + 0(7 3 2 )) . (11) 

We rewrite the temperature T in terms of the deviation from the conductive profile 
(11) neglecting its contribution of 0(jI), 

= T - T cond = T-T -R(-z-^{z 2 - i)) . (12) 

We then obtain as the final dimensionless equations 

1 „ fviW „ , /. , „ 0, 



Pr 



d t v t + Vj dj (-)) = -d tP + 5 l3 (l + 7l (-2z + ^)) © + 



+ a 



+ 2Qe lj3Vj , (13) 



= o, (i4) 

d t e + ^djQ = —djixdje) - l3 d z e -r—(i + -y 3 z). (15) 

p p c p p 

The dimensionless boundary conditions are 

v(x, y, z, t) = 0(x, y, z, t) — at z — ±-. (16) 
The nondimensionalized fluid parameters (7)-(10) read now: 

P(e) =i- 7o (-z + |), (17) 

1/(6) =l + 72 (-^ + |), (18) 
A(0) =1 + ^(^ + 1), (19) 

Cp( e) = l + 74 (-z + |). (20) 

We consider the non-Boussinesq effects to be weak and keep in all material 
properties only the leading-order temperature dependence beyond the Boussinesq 
approximation. Therefore the 71-term appears explicitly in (13), while in all other 
terms it would constitute only a quadratic correction just like the terms omitted in (7)- 
(10). Correspondingly, we expand the denominators in (13,15) that contain material 
properties to leading order in the 7 j. In analogy to [4] (eq.(6.7)), we further omit 
non-Boussinesq terms that contain cubic nonlinearities in the amplitudes Vi or 0, as 
they arise from the expansion of the advection terms Vjdj(vi/p) and (vj/p)djQ, once 
the temperature-dependence of the density is taken into account. Since we will be 
considering Rayleigh numbers up to twice the critical value, which implies enhanced non- 
Boussinesq effects, these approximations may lead to quantitative differences compared 
to the fully non-Boussinesq system, even though the temperature-dependence of the 
material properties themselves may quite well be described by a linear (or quadratic in 
the case of the density) approximation. 
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3. Computational methods 

We solve the NONBE (13,14,15) with top and bottom boundary conditions given by (16) 
and the material parameters by (17,18,19,20) numerically using a number of approaches. 

As usually, the stability properties of the patterns are determined by Galerkin 
methods. We assume an infinite extension of the layer in the lateral directions, which 
is captured by a Fourier expansion on a hexagonal lattice. The Fourier wave vectors 
q are constructed as linear combinations of the hexagonal basis vectors bi = q(l, 0) 
and b 2 = <?(l/2, V3/2) as q = mbi + nb 2 with the integers m and n in the range 
\m\ + | n | < n q . The largest wavenumber is then n q q and the number of Fourier modes 
retained is given by 1 + 6^j=2.7- Typically we used n q = 3. The top and bottom 
boundary conditions are satisfied by using appropriate combinations of trigonometric 
and Chandrasekhar functions in z [17, 4]. In most of the computations we used n z = 6 
modes for each field in eq.(13,14,15). The linear analysis yields the critical Rayleigh 
number R c ('Ji) as well as the critical wavenumber q c ('Ji)- In the large Pr-number 
regime, to be considered in this paper, the bifurcation to convection is stationary. Note 
that q c (li), Rc(li) show 0{p(f) corrections compared to the standard Boussinesq values 
R c = 1708 and q c = 3.118 for % = 0. 

To investigate the nonlinear hexagon solutions, we start with the standard weakly 
nonlinear analysis to determine the coefficients of the three coupled amplitude equations 
for the modes making up the hexagonal pattern. This allows us to get a first insight 
into the Kiippers-Lortz instability of the rolls, and for the weakly non-Boussinesq case 
the transition from hexagons to rolls as well as the transition to oscillating hexagons. 
To obtain the fully nonlinear solutions requires the solution of a set of nonlinear 
algebraic equations for the expansion coefficients with respect to the Galerkin modes. 
This is achieved with a Newton solver for which the weakly nonlinear solutions serve 
as convenient starting solutions. The fully nonlinear solutions are then tested for 
amplitude and modulational instabilities following the standard methods. Alternatively, 
we solve the NONBE by direct numerical simulation to describe in general the temporal 
evolution of the pattern. While this allows for a precise check of the Galerkin results, 
its main purpose is to capture the complex spatio-temporal dynamics resulting from 
the instabilities of the hexagonal pattern, which are not accessible in the Galerkin 
approach. The simulation code involves slight extensions of a previous, well-proven 
version [21, 22, 19] to include the effect of the non-Boussinesq effects. 

4. Results 

Instead of extensive parameter studies, we concentrate in this paper on a specific, 
interesting scenario that should be experimentally realizable. We focus our investigation 
on fluids with a fairly large Prandtl number (water). For smaller Prandtl numbers our 
previous results using a Swift-Hohenberg equation and also some initial test runs using 
the Navier-Stokes equations indicate that rolls [4], which may be nucleated in defects 
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10 


8.728 


0.0012 


0.5888 


-0.2760 


0.0345 


-0.0038 



Table 1. Values for the non-Boussinesq coefficients 7 4 and the Prandtl number 
Pr for water at onset as a function of the mean temperature and the temperature 
difference. 
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12 


2.81 


1708 


1706 


3.12 


3.12 


150.7 





65 


12 


10 


6121 


6080 


5.27 


5.39 


64.6 


0.600 



Table 2. Dependence of the critical Rayleigh number and of the critical wavenumber 
on the rotation rate for a layer thickness of d = 0.46cm. The critical Rayleigh numbers 
in the Boussinesq approximation = 0) and for 7$ = 7^(AT = AT C ) are denoted by 
R c o and R CJ , respectively. Similarly for the critical wavenumber q c0n . The physical 
rotation rate for this cell thickness is given by u> and the vertical thermal diffusion time 
by r th . 



or at the wall, take over the hexagon pattern already for relatively small values of 
the Rayleigh number. Furthermore, we have confined ourselves in this work on one 
representative rotation rate (Q = 65, see below). 

To obtain strong non-Boussinesq effects the layer would have to be taken quite thin. 
However, then the dimensionless rotation rates obtained for a given physical rotation 
rate are relatively low. Therefore we decided to consider water near its density maximum 
as the convecting fluid. Strong non-Boussinesq effects at larger layer thickness are also 
obtained with glycerin. However, since the rotation rate is made dimensionless using 
the viscous diffusion time d 2 /u the extremely large viscosity v of glycerin yields only 
small dimensionless rotation rates. The second row in Table 1 gives the non-Boussinesq 
parameters used in all computations in this paper. For Q = 65 they correspond to a 
critical temperature difference of 10°C across the layer at a mean temperature of 12°C 
Note that for a given thickness of the fluid layer the critical Rayleigh number increases 
significantly with the rotation rate. Therefore, the coefficients 7^, which give the non- 
Boussinesq effects at onset, also increase strongly with the rotation rate. Consequently, 
while in the absence of rotation the system is essentially Boussinesq for the chosen layer 
thickness of d = 0.46cm (first line in table 1) it is strongly non-Boussinesq for = 65. 

We first determine numerically the linear stability of the stationary hexagon 
patterns as a function of the wavenumber and the control parameter e = (R — 
Rc(li)) / Rcili) for f2 = 65 and the corresponding non-Boussinesq coefficients 7^ as given 
in table 1. As indicated before, the 7^ are linear in the temperature difference AT and 
are therefore given by 7^ = 7^ (1 + e). All results were obtained with n q = 3 and n z = 6. 
In the parameter regime investigated these cut-off parameters are sufficient. As shown 
in Fig.l the hexagon pattern can undergo a variety of linear instabilities in the nonlinear 
regime. Close to threshold the dominant instability is a long-wave instability, i.e. its 
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growth rate vanishes for vanishing modulation wave number f3 (Floquet parameter) but 
grows quadratically with (3 for small (3. Very close to onset it can be steady (black 
circles). For slightly larger e, however, it becomes oscillatory (red squares). As the 
control parameter is increased a steady instability at finite Floquet parameter limits the 
stability region of the hexagons for large wavenumbers (blue triangles). Increasing the 
control parameter, on the low-wavenumber side the stability region becomes limited by 
an oscillatory instability that sets in with a finite Floquet parameter (green diamonds). 
For non-Boussinesq coefficients that differ slightly from those of table 1 (by (9(20%)) an 
oscillatory instability with vanishing Floquet parameter determines the stability limit 
for low q and e ~ 1. For the parameters in table 1 it is, however, preempted by the 
Hopf bifurcation at finite modulation wavenumber. 

The spatially homogeneous and the short-wave oscillatory instabilities are similar 
to the oscillatory instability of hexagons that is found within the framework of three 
coupled Ginzburg-Landau equations [9, 10, 5, 6]. There it replaces the steady instability 
of hexagons that in the absence of rotation leads to an unstable mixed-mode solution. 
With rotation it turns into a supercritical Hopf bifurcation at which the oscillating- 
hexagon solution branches off the stationary hexagons. 




4.4 4.6 4.8 5 5.2 5.4 5.6 



Wavenumber q 

Figure 1. Stability limits of stationary hexagon pattern for the parameters given in 
Table 1 (n q — 3, n z — 6, except for dark green diamond, which is for n q = 5, n z — 8): 
long- wave steady (black circle), long- wave oscillatory (red squares), oscillatory short- 
wave instability (green diamonds), steady short-wave instability (blue triangles). 

Numerical simulations support the expectation that the modes corresponding to 
the Hopf bifurcation are similar to the oscillating hexagons. In a small computational 
domain of length L = 2 ■ 2tc /q, which contains only four convection cells, the oscillatory 
side-band instability is suppressed and the homogeneous Hopf bifurcation arises when 
the parameters are chosen sufficiently far beyond the stability limit corresponding to 
the short-wave oscillatory instability (green diamonds in Fig.l). This instability leads 
to an elliptic deformation of the convection cells which oscillates in time, giving the cells 
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the appearance as if they were rotating or in fact whirling. In sufficiently large domains 
the side-band instability leads to a spatial modulation of these oscillations as shown in 
Fig.2a (L — 8 • 2n/q with 64 x 64 Fourier modes). 




Figure 2. Snapshots of the pattern for e = 1.0, L = 8 x 2n/q with q = 4.5 (see also 
movie sn76.T.mpg). Initial growth of modulational instability at t = 7.5 (left panel) 
and t = 7.9 (middle panel). Disordered state during an intermediate phase at t = 51.3 
(right panel). Black denotes hotter rising fluid. 



As the instability develops and the oscillation amplitude grows the regular spatial 
structure is destroyed and defects arise in the pattern as pictured in Fig.2c. This 
temporal evolution is illustrated in the movie sn76 . T . mpg. While the apparently chaotic 
dynamics appear to be reaching a statistically steady state, they do, however, not persist 
and around t = 130 a relative rapid transition to an ordered, stationary hexagon pattern 
occurs. It has a slightly larger wave number q, which is in the stable regime (cf. Fig. 
1). During this chaotic transient the overall disordered pattern intermittently exhibits 
relatively large domains of ordered hexagonal patterns. The duration of the transient 
increases with e as illustrated in Fig. 3, reflecting presumably the narrowing of the band 
of stable wavenumbers with increasing e. To monitor the activity of the pattern we 
consider as a representative measure the temperature field in the mid plane of the layer 
and plot the 'whirling' activity W given by the Li-norm of the time derivative of its 
Fourier modes $ normalized by the Li-norm of the Fourier modes themselves. As e is 
increased the fluctuations in this quantity decrease, which is expected to reduce the 
probability of excursions to one of the linearly stable stationary states, and eventually 
for e = 1.1 our simulations did not end up in a stationary pattern for times as large as 
t = 400. 

The duration of the transients is presumably related to the size of the system. 
We therefore performed also simulations in larger systems of size 16 • 2-n/q (128 x 128 
Fourier modes). For this system size the dynamic state persists for quite a long time 
{tmax ~ 300r t ft « 5hrs) even for control parameters as low as e = 0.87 (compare with 
the transients for the smaller system shown in Fig. 3). This increased duration of the 

| Since the temperature field is dominated by the linear eigenvector with the maximal growth rate this 
allows to reduce the computational overhead. 
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100 200 300 JOO 

Tim.*t 



Figure 3. Dependence of the duration of the transient on e for L = 8 x 2ir/q with 
q — 4.5 as quantified by the whirling activity W (see text). In physical units t = 400 
corresponds to about 7 hours. 



transient is due to the longer time it takes for the defects in the pattern to annihilate. For 
e = 1.0, a snapshot of the pattern and of the corresponding temporal derivative is shown 
in Fig.4a,b. The latter makes the spatially intermittent behavior of the dynamics more 
apparent. The dynamics are illustrated in the two corresponding movies (sn79h.T.mpg, 
sn79h . D . mpg) . 
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Figure 4. Disordered state for e = 1.0 and L = 16 x 27r/g with g = 4.5. Snapshot 
of the pattern (blue denotes hotter rising fluid); see also movie sn79h.T.mpg (left 
panel) . Middle panel shows the difference between two successive snapshots (separated 
by At — 0.1) corresponding to the pattern shown in left panel (see also movie 
sn79h . D . mpg) . Blue (red) denotes locations with growing (decreasing) temperature. 
Right panel: identification of the minima with lattice nodes. Blue (red) circles have 5 
(7) nearest neighbors. Up- (down-) triangles have 4 (8) neighbors. 

For the characterization of weakly disordered hexagon patterns it is often 
advantageous to locate defects in the pattern by demodulating the pattern using 
the (three) main wavevectors comprising the pattern (e.g. [23, 7]). This allows the 
identification of the penta-hepta defects that are characteristic of weakly disordered 
hexagon patterns and which consist of correlated convection cells with five and seven 
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neighbors, respectively. Obviously, this approach requires that the Fourier spectrum of 
the pattern exhibits well-defined peaks at the three hexagon modes. 

In the strongly disordered patterns obtained in our simulations the spectrum does 
not exhibit clear peaks. Demodulation is therefore not possible and other techniques 
for the analysis of the patterns have to be brought to bear. We use a triangulation 
of the pattern to locate the defects. To that end we map the convection pattern to a 
point lattice by identifying the local minima in the temperature with nodes in a lattice. 
To avoid spurious minima we require that the minima fall below a certain threshold. 
The coordination number of each node is then obtained using a Delaunay triangulation, 
which is the unique triangulation obtained by the requirement that the circle through 
the corners of each triangle do not contain any other node of the pattern. We use 
the program triangle [24] to perform the triangulation. Figure 4c shows the lattice 
of the snapshot shown in Fig. 4a with the local coordination numbers obtained from 
the triangulation. Black diamonds indicate nodes with six neighbors, i.e. the center 
of a hexagonal cell. Nodes with seven (i.e. heptagons) and five neighbors (pentagons) 
are marked with red and blue circles, respectively, while up- and down-triangles denote 
nodes with four and eight neighbors, respectively. For less disordered patterns pentagons 
and heptagons are often strongly correlated and can be associated with each other in 
pairs. These pairs correspond then to penta-hepta defects and can be identified with 
pairs of topological dislocation defects in two of the three modes that make up the 
hexagonal pattern (see e.g. Fig.25 in [19]). Isolated pentagons or heptagons, however, 
do not correspond to single dislocations. 

In previous analyses of defect-dominated spatio-temporal chaos it was found that 
the distribution function of the number of defects in the patterns gives some insight 
into the dynamics of the system. More specifically, in the defect chaos obtained in the 
complex Ginzburg-Landau equation [12], in electroconvection of nematic liquid crystals 
[14], and in undulation chaos in inclined-layer convection [15, 16] the statistics were 
found to be close to a Poisson-type distribution. Such a distribution arises if defects of 
opposite topological charge are created randomly in pairs with a fixed probability and 
are annihilated with a rate that is proportional to the square of the density of the defects. 
The Poisson statistics therefore suggest that the defects move like uncorrelated random 
walkers. In contrast, for the penta-hepta defect chaos that was found in hexagon patterns 
in the presence of rotation within the framework of a Swift-Hohenberg-type equation 
the distribution function was found to be considerably broader [7]. The origin of this 
broadening was identified to be the induced nucleation of dislocations by penta-hepta 
defects through which also the creation rate depends on the defect density [7, 25]. 

Using the triangulation procedure, we have measured the distribution functions for 
the number of pentagons and heptagons. For the state corresponding to Fig.4 they are 
presented in Fig.5a,b, where the blue circles and red triangles denote the distribution 
functions for pentagons and heptagons, respectively. As mentioned previously, in 
strongly disordered patterns pentagons and hexagons are not necessarily bound in penta- 
hepta defects. Correspondingly, the number of pentagons does not always agree with 
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Figure 5. a) Probability distribution function for the number of pentagons (blue) and 
heptagons (red) for e — 1 and corresponding best fit to the Poisson distribution (21) 
(dashed lines). 

b) Probability distribution function for the number of rectangles (blue), octagons (red), 
and for the difference (V n5 -n 7 ) between the number of pentagons and heptagons (solid 
line). Parameters as in Fig. 4. 



the number of heptagons. The distribution function for their difference is given by the 
dashed line in Figure 5b along with the distribution function for rectangles (blue circles) 
and octagons (red triangles). 

As a first step it appears reasonable to assume that the dynamics of the pentagons 
and heptagons are captured by a kinetic model that is analogous to the description of 
dislocations [12]. For periodic boundary conditions [15] this leads then to the standard 
Poisson distribution for each species, 

nN) = /(2^)r(i + Ar (21) 

where a =< N 2 > and I(2y/a) is the modified Bessel function. In some simulations we 
observe that a given whirling cell creates a new convection cell nearby, which implies 
the creation of defects. The newly created cell then often merges with the whirling 
cell. Conversely, defects in the pattern induce whirling activity. This mutual nonlinear 
enhancement of defects and whirling activity suggests that the existence of defects may 
lead to the creation of further defects, somewhat similar to the induced nucleation in 
penta-hepta defect chaos [7]. As shown by the dashed lines in Fig.5a, for the state 
of Fig.4 the Poisson distribution (21) does not fit the pentagon and the heptagon 
distributions very well. The distributions therefore suggest that the creation rates 
for the pentagons and heptagons are not independent of their respective densities. In 
some preliminary parameter studies involving non-Boussinesq coefficients that may not 
correspond to a realistic convection cell with water as fluid, we have found distributions 
that are fitted quite well by the Poisson distribution and others that are even broader 
than the distributions obtained in Fig.5. This suggests that a mechanism akin to 
induced nucleation may be operative in this system, with its significance depending 
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on the specific parameter values [26]. A study of the dependence of the distribution on 
the physical parameters of the fluid is beyond the scope of the present paper. 

Fig.4b shows that in the disordered state not all cells are active and at times sizable 
patches of relatively ordered quiescent domains arise. To characterize this intermittent 
behavior we consider the distribution function for the time derivative of the temperature 
as obtained from all snapshots of the type shown in Fig. 4b. The logarithmic scale of Fig.6 
reveals an exponential decay of the distribution function for e = 1.0. This exponential 
dependence is consistent with the visually apparent spatio-temporal intermittency of 
the activity of the pattern. The pronounced peak exhibited for e = 0.87 represents the 
large spatial domains in which the hexagons are ordered and quiescent. 




Difference between successive snapshots 



Figure 6. Probability distribution function for the difference between successive 
snapshots (At = 0.1) for L — 16 x 2-k/8 and q = 4.5. 

In all simulations so far we have employed periodic boundary conditions. In 
experiments boundary effects may strongly perturb the patterns. We therefore have also 
performed some simulations in which we mimic a circular container by applying a strong 
radial subcritical ramp in the Rayleigh number that suppresses any convection outside 
a certain radius. In previous simulations [19] this was seen to provide a reasonable idea 
of the impact of boundaries. Experimentally, it has been found that boundaries tend 
to induce defects and may also promote the formation of rolls, which could invade the 
system and replace the hexagonal pattern [27]. For the parameter regimes investigated 
here we find that no patches of rolls invading from the side replace the hexagons. Instead 
the precession of the cells near the boundary and the associated formation of defects 
triggers persistent dynamics in the whole convection cell even for parameter values for 
which in the absence of the boundaries the pattern eventually becomes ordered and 
stationary. Fig. 7a shows a snapshot for e = 0.87 with the associated temporal evolution 
displayed in the movie sn74 . c . T . mpg. A snapshot for a larger system size and e = 1.0 is 
shown in fig. 7b (movie sn79h . c . T . mpg) . We therefore expect that the defect-dominated 
spatio-temporal chaos will also be accessible in experiments. 

We have not investigated in detail the nonlinear evolution ensuing from the steady 
instability at finite wavenumber that limits the wavenumber band on the high-g side. 
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Figure 7. Snapshot of defect-dominated pattern in the presence of a radial ramp 
mimicking a circular container, a) e = 0.87, L = 8 • 2ir/q with q = 4.5. b) e = 1.0, 
L = 16 • 27r/g with q = 4.5. Temporal evolution shown in movies sn74 . c . T . mpg and 
sn79h.c.T.mpg.) 

In preliminary simulations the instability led to the formation of defects in the pattern, 
which eventually, however, disappeared yielding a stationary, regular hexagonal pattern. 

5. Conclusion 

In this paper we have considered hexagonal non-Boussinesq convection focusing on 
rotating systems, in which the chiral symmetry is broken. The main motivation for 
this work stems from the interesting chaotic states that had been identified in previous 
investigations of the stability of hexagonal patterns and of the dynamics resulting from 
their instabilities. The chaotic states were associated with the secondary instability to 
oscillating hexagons [5, 6] and with the induced nucleation of dislocations [23, 7, 25], 
respectively. These results were obtained within order-parameter and Ginzburg-Landau 
models. To make closer contact with experimental systems we have investigated in this 
paper the effect of rotation on realistic hexagonal convection in a thin layer of water in 
a temperature regime near the anomalous density maximum. 

Using a numerical stability analysis and direct numerical simulations of the Navier- 
Stokes equations we have identified an oscillatory side-band instability, which leads 
to spatially modulated whirling hexagons that in turn become disordered in space and 
exhibit irregular persistent dynamics, which coexists with linearly stable regular hexagon 
patterns. 

Coexistence of spatio-temporal chaos with ordered patterns has been observed 
previously in roll convection at low Prandtl numbers [1, 28]. There the spiral defect 
chaos is sustained by the mean flow, which compresses the rolls so as to push them 
locally across their stability limit with respect to the skew-varicose instability, which 
then leads to the formation of defects [29]. This mechanism for sustaining dynamics is 
reminiscent of that observed earlier in a much smaller systems [30, 31]. 
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The method employed in [29] to analyze the disordered pattern, which extracts the 
local wavenumber and orientation of the rolls, is not applicable to the strongly disordered 
hexagonal patterns found in our simulations. We therefore use Delaunay triangulation 
to measure the distribution functions for the number of pentagons and heptagons. Both 
show deviations from a squared Poisson distribution, which suggests that their dynamics 
(including their creation and annihilation) exhibit correlations that are absent in the 
dynamics of dislocations in the complex Ginzburg-Landau equation [12] or in inclined- 
layer convection [15, 16]. It is, however, not quite clear whether a kinetic description of 
the pentagons and heptagons similar to that of dislocations or of penta-hepta defects 
[12, 15, 7] is in fact appropriate. To characterize the intermittent appearance of ordered 
domains in the chaotic states we have determined the distribution function for the 
temporal derivative of the temperature field. It exhibits an exponential decay reflecting 
the spatio-temporal intermittency. For smaller Rayleigh number a pronounced peak 
appears, which reflects the large ordered, quiescent domains. 

In this paper we have focused on one set of fluid parameters (cf. table 1). For 
these values the oscillatory instability with finite modulation wavenumber preempts a 
spatially homogeneous oscillatory instability. We found, however, that the situation can 
be reversed if the non-Boussinesq parameters are changed only slightly. In these cases, 
which may not correspond to a realistic convection cell with water as fluid, we have 
found defect-chaotic states with distribution functions that are significantly wider than 
the squared Poisson distribution suggesting again a nonlinear interaction between the 
whirling of the hexagons and the formation of defects, as well as distribution functions 
that are well described by the squared Poisson distribution. In addition, we have 
identified a state in which the whirling is spatially and temporally intermittent in the 
form of bursts while the underlying hexagonal planform has essentially no defects [26]. 
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